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1.  Introduction 

MicroElectroMechanicalSystems  (MEMS)  are  often  characterized  by  structures  that  are  a  few 
microns  in  size,  separated  by  micron-sized  gaps.  At  these  sizes  and  gaps,  viscous  air  damping 
dominates  other  energy  dissipation  mechanisms  at  atmospheric  pressure,  thus,  affecting  the 
dynamic  behavior  of  the  devices.  Also,  the  number  of  such  gaps  can  be  very  large  in  a  MEMS 
device.  Hence,  fast  and  accurate  behavioral  simulation  of  the  dynamics  of  MEMS  systems  necessi¬ 
tates  accurate,  low-order  damping  models.  Further,  to  be  useful  in  the  design  flow,  these  models 
need  to  be  parameterized,  so  that  they  can  be  used  to  model  damping  for  a  wide  range  of  structure 
size,  gap,  frequency  and  amplitude  of  motion. 


Viscous  damping  due  to  the  fluid  surrounding  the  moving  structures  can  be  of  two  types: 
squeeze  film  damping  and  lateral  damping.  In  squeeze  film  damping,  there  are  two  parallel  struc¬ 
tures  as  shown  in  Fig.  1.1.  They  are  in  relative  motion  perpendicular  to  the  plane,  thus  squeezing 
the  air  film  between  the  structures.  The  air  counters  the  motion  by  exerting  an  opposing  force.  Fig. 
1 .2  shows  the  case  of  lateral  damping.  Here  the  motion  of  the  plate  is  parallel  to  the  plane  of  the 
structure.  The  resulting  fluid  flow  causes  energy  dissipation  and  exerts  an  opposing  force  which  is 
a  function  of  the  velocity  of  the  moving  plate,  the  dimensions  of  the  plate,  its  distance  from  other 
parallel  structures,  and  the  viscosity  of  the  fluid. 
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Fig.  1.1:  Cross  section  of  a  squeeze  film  between  a  moving  plate  and  the  substrate 


This  report  presents  lumped  parameter  models  for  squeeze  film  and  lateral  damping.  The 
accuracy  of  these  models  down  to  MEMS  scaled  structures,  and  their  parameterized  nature  makes 
them  useful  for  design  and  simulation  of  MEMS  devices.  With  the  edge-effect  contribution  to  the 
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Fig.  1.2:  Cross  section  of  a  moving  plate  that  is  laterally  damped 


damping  force  becoming  significant  in  micron-scaled  structures,  emphasis  is  laid  on  including 
these  effects. 

In  chapter  2,  squeeze  film  damping  model  is  discussed.  Finite  element  simulations  of 
squeeze  film  are  presented.  The  derivation  of  the  behavioral  model  for  squeeze  film  from  existing 
analytical  model  [1]  and  its  extension  to  include  the  edge  and  the  finite-size  effects  are  discussed. 
The  significance  of  the  contribution  of  these  effects  for  MEMS  scaled  devices  is  established.  The 
edge-effects  enhanced  behavioral  model  is  evaluated,  for  varying  frequency,  amplitude  of  plate 
motion  and  the  aspect  ratio  of  the  plate  by  comparing  it  with  FEM  simulations.  Simulation  results 
for  a  square  wave  input  are  presented. 

In  chapter  3,  modeling  of  the  fluid  flow  invoked  by  a  laterally  moving  structure  and  the  deri¬ 
vation  of  the  behavioral  model  for  lateral  damping  are  discussed.  The  accuracy  of  the  behavioral 
model  is  quantified  by  comparison  with  the  finite  element  simulations.  Simulation  results  for  a 
square  wave  input  and  for  varying  aspect  ratio  and  oscillation  frequency  are  discussed. 

Chapter  4  establishes  the  use  of  the  presented  damping  models  in  system  level  simulation  and 
design,  by  use  in  the  nodal  simulation  of  a  bandpass  filter  and  a  lateral  resonator.  The  simulation 
results  are  compared  with  the  experimental  results  which  establishes  the  accuracy  of  the  damping 
models. 


Finally,  chapter  5  concludes  the  overall  work  and  discusses  some  future  work  that  could  be 
done  in  this  direction. 
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2.  Squeeze  film  damping 


2.1  Introduction 

When  a  fluid  film  between  two  parallel  surfaces  is  squeezed  by  the  relative  motion  of  the  sur¬ 
faces  towards  each  other,  a  force  that  opposes  the  motion  of  the  structures  is  produced.  This  type 
of  energy  dissipation  mechanism,  called  squeeze  film  damping,  is  encountered  very  commonly  in 
MEMS  devices.  A  squeeze  film  is  characterized  by  physical  dimensions  like  width,  length  and 
thickness  of  the  fluid  film  that  is  squeezed.  Other  parameters  that  determine  the  damping  are  the 
ambient  pressure,  the  temperature  of  device  operation,  viscosity  coefficient  of  the  fluid  and  the 
velocity  of  the  oscillating  plate. 

Example  applications  include  MEMS  systems  that  have  a  proof  mass  oscillating  vertically 
above  a  static  bottom  plane  or  a  differential  comb  drive.  Damping  of  the  proof  mass  in  a  z-axis 
accelerometer  is  a  fundamental  consideration  in  the  design  of  the  accelerometer.  In  differential 
comb  drives,  the  squeezed  air  between  the  stator  and  the  rotor  fingers  contributes  most  of  the 
damping  in  the  differential  comb  drives.  Fig.  2.1  shows  the  squeeze  film  between  the  stator  and 
rotor  fingers  of  a  differential  comb  drive.  The  stator  finger  moves  in  the  direction  marked  in  Fig. 
2.1  (b).  The  thickness  of  the  structure  is  equal  to  the  width  of  the  squeeze  film  (typically  a  few 
microns);  the  overlap  distance  of  the  stator  and  rotor  fingers  is  equal  to  the  length  of  the  squeeze 
film  (tens  of  microns);  and  the  equilibrium  gap  between  the  stator  and  the  rotor  fingers  is  equal  to 
the  gap  between  the  top  and  the  bottom  plate  (typically  a  couple  of  microns). 


Rotor. 


Rotor 


Displacement 
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Stator 
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(a)  (b) 

Fig.  2. 1  (a):  SEM  picture  showing  the  stator  and  the  rotor  fingers  of  a  differential  comb  drive 
(b):  Picture  showing  the  squeeze  film  between  the  stator  and  rotor  fingers 
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2.2  Squeeze  film  model 
2.2.1  Squeeze  film  physics 


The  gas  pressure  distribution  across  the  moving  plate,  as  a  result  of  squeeze  film  effect,  is 
governed  by  the  Reynolds  equation.  Under  isothermal  conditions,  the  pressure  (that  is  close  to 
ambient  pressure)  is  inversely  proportional  to  the  volume  occupied  by  the  gas,  or,  directly  propor¬ 
tional  to  the  density  of  the  gas.  When  the  plates  are  parallel  and  the  relative  displacement  of  the 
moving  plate  is  small  as  compared  to  the  distance  between  the  plates,  the  second-order,  non-linear 
partial  differential  equation  simplifies  to  the  linearized  Reynolds  equation  (EQ.  2.1)  and  to  EQ.  2.2 
under  isothermal  conditions  [2,  3]. 
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where  p  is  the  density  of  the  gas,  p  is  the  change  in  pressure  from  the  ambient  pressure,  Pa  is  the 
ambient  pressure,  g  is  the  gap  between  the  plates,  r)  is  the  viscosity  coefficient  of  air  and  z  is  the 
vertical  displacement  of  the  moving  plate,  x  and  y  are  the  cartesian  coordinates. 


2.2.2  Gas  rarefaction  effects 

The  squeeze  film  sizes  that  are  found  in  MEMS  devices  may  result  in  a  condition  where  the 
fluid  film  can  no  longer  be  considered  to  be  within  continuum  limits.  As  the  film  thickness  is 
reduced,  the  mean-free  path  of  the  molecules  becomes  significant  as  compared  to  the  critical 
dimension,  which  is  the  squeeze  film  thickness  in  this  case.  This  regime  of  operation  can  be  given 
by  a  dimensionless  quantity  called  the  Knudsen  number,  Kn,  which  is  the  ratio  of  the  mean  free 
path  of  the  gas  molecules  under  consideration  X,  to  the  gap  between  the  plates,  g.  When 
Kn  <  1/100,  the  fluid  can  be  considered  to  be  in  the  continuum  limits  and  the  no-slip  assumption 
that  the  plate  velocity  and  the  fluid  velocity  at  plate  surface  is  equal  remains  valid.  On  the  other 
end,  if  Kn>  1 ,  then  molecular  regime  of  operation  should  be  considered.  Slip  conditions  begin  to 
occur  when  X  becomes  a  fraction  of  g  (>1/100).  The  mean  free  path  of  air  is  about  70  nm  at  stan¬ 
dard  temperature  and  pressure.  The  gaps  that  are  commonly  encountered  in  squeeze  film  damping 
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can  be  a  couple  of  microns.  So  slip  conditions  have  to  be  taken  into  account.  At  lower  pressures  Kn 
is  higher  because  X  increases,  and  the  region  of  operation  shifts  more  and  more  towards  the  molec¬ 
ular  regime. 

The  effect  of  slip-conditions  can  be  modeled  by  using  a  value  of  effective  viscosity  which  is 
lower  than  the  actual  viscosity  of  the  fluid.  There  are  many  existing  equations  to  model  the  viscos¬ 
ity  of  the  fluid  flowing  in  a  narrow  channel  as  summarized  in  [3].  We  use  the  value  pointed  to  by 
Veijola  et  al.  [3].  The  equation  for  effective  viscosity  is  given  by: 

T]  EQ.  2.3 

^  1  +  9. 638K[ 159 

where  r^is  the  effective  viscosity,  r|  is  the  fluid  viscosity,  Kn  is  the  Knudsen  number.  Fig.  2.2 
plots  the  effective  viscosity  as  a  function  of  gap  between  the  plates,  the  decrease  in  the  effective 
viscosity  can  be  seen  with  decreasing  gap. 


Fig.  2.2:  Effective  viscosity  as  a  function  of  gap 

2.2.3  Squeeze  force:  Solution  to  compressible,  linearized  Reynolds 
equation 

The  linearized  Reynolds  equation  (EQ.  2.2)  can  be  solved  under  compressible  conditions  for 
sinusoidal  oscillation  of  the  a  rectangular  plate.  The  solution  consists  of  a  spring  force  in  phase 
with  the  plate  displacement  and  a  damping  force  in  phase  with  the  plate  velocity[Blechs],  The  inte- 


5 


gral  of  the  gauge  pressure  over  the  area  of  the  plate  gives  the  following  relations  for  the  damping 


and  the  spring  forces  [1]: 


Fd  _  64oPawl  ^  m2+n2c 2 

z  71  9  ~(mnf{[nf+rfc*]*+<fA4} 


EQ.  2.4 


Fs  64<52Pawl  ^  1  EQ.  2.5 

n*g  mnftfnf  +  nzc2  ]2  W M4  } 

odd 


where  F^and  Fs  are  the  damping  and  spring  forces,  /  is  the  length,  w  is  the  width  and  c  is  the  aspect 
ratio  of  the  plate  (width/length),  z  is  the  amplitude  of  motion,  g  the  gap  and  a  is  the  squeeze  num¬ 
ber  given  by: 

(!) 

EQ.  2.6 


o  = 


12r\effw  (D 


where  co  is  the  angular  frequency  of  oscillation. 


2.3  Finite  element  simulation  of  squeeze  film 

Finite  element  simulations  of  the  dynamics  of  squeeze  films  are  performed  using  a  code 
which  solves  three-dimensional  Navier-Stokes  equation.  One  such  commercial  code  used  in  this 
thesis  is  CFD-ACE  [4,  5].  Transient  simulations  of  a  plate  oscillating  vertically  above  an  infinite 
substrate  are  performed.  The  damping  force  on  the  oscillating  plate  is  determined  by  the  integra¬ 
tion  of  the  gauge  pressure  (pressure  deviation  from  the  ambient  pressure)  over  the  surface  of  the 
oscillating  plate. 


To  verify  the  finite  element  simulations,  a  well  known  and  analyzed  problem  was  solved.  A 
square  plate  of  size  1mm  on  a  side,  oscillating  with  an  amplitude  of  0.1  pm  and  a  frequency  of  1 
kHz  was  chosen.  Fig  2.3  shows  the  damping  force  with  time  for  a  gap  of  2  pm.  The  damping  coef¬ 
ficient  can  be  extracted  from  the  peak  force.  The  analytical  equation  for  the  damping  coefficient 
for  a  square  plate  of  length  l  on  a  side  and  a  gap  g,  is  given  by  EQ.  2.7  [6]. 

B  =  0422 —il—  EQ.  2.7 


where  r)eff  is  the  effective  viscosity  coefficient  of  the  fluid. 
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The  values  of  damping  coefficient  obtained  from  the  FEA  method  with  ambient  pressure  set 
at  the  plate  edges  and  the  analytical  method  are  0.96  and  1.0  respectively.  They  match  to  within 
4%  which  verifies  and  quantifies  the  accuracy  of  the  finite  element  simulations. 


Fig.  2.3:  Damping  force  on  a  sinusoidally  oscillating  plate 


The  accuracy  of  the  finite  element  simulations  also  depends  on  the  boundary  conditions  used. 
Two  different  boundary  conditions  can  be  used  in  the  simulation  of  the  squeeze  film  damping:  triv¬ 
ial  boundary  conditions  and  non-trivial  boundary  conditions.  The  advantages  of  each  approach  is 
discussed  next. 


(a):  Squeeze  film  with  TBC  (b):  Pressure  distribution  across  top  plate 

Fig.  2.4:  20  pm  square  plate  oscillating  with  0.1  pm  amplitude  at  1  kHz 
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2.3.1  Trivial  boundary  conditions 

In  trivial  boundary  conditions,  the  squeeze  film  space  is  represented  by  a  3-D  grid  between  an 
oscillating  plate  and  substrate  as  shown  in  Fig.  2.4  (a).  To  exploit  the  symmetry  conditions  for 
reducing  the  number  of  grid  cells  and  the  computer  resources  in  terms  of  memory  and  CPU  time, 
only  a  quarter  section  of  the  squeeze  film  is  used  in  simulations.  Ambient  pressure  boundary  con¬ 
dition  is  used  at  the  plate  edges.  Fig.  2.4.b  shows  the  pressure  distribution  across  a  the  quarter  plate 
when  a  square  plate  of  size  20  pm  is  oscillating  sinusoidally  with  an  amplitude  of  0.1pm  and  a  fre¬ 
quency  of  1  kHz.  It  can  be  noticed  that  the  gauge  pressure  is  forced  to  zero  at  the  plate-edges. 


Ambient  Pressure 


Fig.  2.5  (b):  Pressure  distribution  across  and  around  the  top  plate  (20  pm  X  20  pm) 
oscillating  with  an  amplitude  of  0.1  pm  at  1  kHz 
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2.3.2  Non-trivial  boundary  conditions 

The  use  of  non-trivial  boundary  conditions  are  more  realistic  boundary  conditions  by  allow¬ 
ing  the  solver  to  arrive  at  the  gauge  pressure  at  the  plate  edges.  This  is  done  by  extending  the  con¬ 
trol  volume  around  the  oscillating  plate  and  setting  the  pressure  to  ambient  pressure  at  the 
extended  control  volume  boundary.  Fig.  2.5  (a)  shows  the  system  with  non-trivial  boundary  condi¬ 
tions  with  blocks  of  air  surrounding  the  plate  oscillating  above  an  infinite  substrate.  Fig.  2.5  (b) 
shows  the  pressure  distribution  across  the  quarter  plate.  Notice  that  the  pressure  contours  extend 
beyond  the  plate  edges. 


2.3.3  Trivial  boundary  conditions  vs.  non-trivial  boundary  conditions 

In  practice,  the  pressure  at  the  plate  edges  is  non-zero.  When  trivial  boundary  conditions 
(TBC)  are  used,  this  value  is  forced  to  zero  and  causes  underestimation  of  the  gauge  pressure  dis¬ 
tribution  across  the  plate.  This  results  in  the  underestimation  of  the  damping  force  due  to  squeeze 
action.  This  problem  can  be  mitigated  by  the  use  of  non-trivial  boundary  conditions  (NTBC) 
which  models  the  boundary  conditions  more  realistically. 

The  following  table  provides  the  squeeze  force  predicted  by  the  use  of  NTBC  and  TBC  when 
a  plate  of  varying  size  oscillates  with  an  amplitude  of  0.1  pm  and  frequency  of  1  kHz,  2  pm  above 
the  substrate.  It  can  be  observed  from  the  table  that  the  difference  in  force  predicted  by  TBC  and 
NTBC  decreases  as  the  plate  size  increases.  This  can  be  intuitively  explained  because  the  edge- 
effects  and  the  finite-size  effects  decrease  with  increasing  plate  size  reaching  zero  for  an  infinitely 
large  plate.  The  dominance  of  the  edge  and  the  finite-size  effects  for  sub-100  pm  plates  has 
already  been  experimentally  verified  in  case  of  lateral  viscous  damping.  For  developing  behavioral 


Plate  Size 

TBC 

NTBC 

Difference 

20pm 

109.1  pN 

151.0  pN 

27.8% 

100pm 

62.5  nN 

67.1  nN 

7.6% 

1mm 

0.65  mN 

0.67  mN 

2.9% 

models  for  damping  for  use  in  MEMS-scale  structures,  all  finite  element  simulations  are  done 
using  non-trivial  boundary  conditions  to  take  care  of  the  edge  and  the  finite-size  effects. 
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The  use  of  NTBC  results  in  the  increase  in  the  number  of  grid  cells  from  500  to  5000  as  a 
result  of  which,  a  typical  transient  simulation  time  increases  from  5  minutes  in  TBC  to  60  minutes 
in  NTBC  on  a  360  MHz  Sun  Ultra-Sparc.  The  memory  requirements  of  the  simulation  is  5  MB  in 
TBC  and  45  MB  in  NTBC.  Though  there  is  an  increase  in  the  use  of  computer  resources  in  terms 
of  CPU-time  and  memory,  this  provides  the  solver  with  more  realistic  boundary  conditions  to  give 
accurate  results. 

2.4  Lumped  parameter  model 

A  lumped  parameter  model  consisting  of  parallel  branches  of  series  connected  spring  and 

damper  elements  can  be  derived  to  model  the  damping  force  given  by  EQ.  2.4  and  EQ.  2.5,  with 

force  as  the  through  variable  and  velocity  as  the  across  variable.  To  derive  the  values  of  the  spring 

and  the  damper  elements,  we  consider  the  following  spring-damper  element  ladder: 

Force  F 

— ►  rrrvm _ q= _ 

+  velocity  v 

Fig  2.6  A  branch  with  series  connected  spring  and  damper  elements 


The  velocity  across  the  terminals  of  the  ladder  is  the  sum  of  the  relative  velocities  across  the 
spring  and  the  damper  elements  (EQ.  2.7  and  EQ.  2.8).  Considering  the  through  variable  and  the 
across  variable  to  be  vectors  of  angular  frequency  to,  the  final  relation  is  described  in  EQ.  2.11. 


v  =  v  . 

spring  damper 


-  F  d 
v  =  — + — i 
B  dt 


. -  F  .(P^ 


2  ,  ,2 


kB  co 


k'  +q)2B2 


+  / 


Bk  co 


EQ.  2.8 
EQ.  2.9 

EQ.  2.10 
EQ.  2.11 


k  +0)  B 


z  o2 


The  real  part,  in  phase  with  the  displacement  vector,  is  the  spring  component  and  the  imagi¬ 
nary  part  is  the  damping  component.  Upon  comparing  the  imaginary  and  real  parts  of  EQ.  2.11 


to 


with  the  force  terms  in  EQ.  2.4  and  EQ.  2.5  respectively,  the  values  of  the  spring  and  the  damper 


elements  is  obtained  for  a  given  pair  of  odd  integers  m,  n  [7]. 
6!wPa 

(mnf  n4  (g+z) 

B  76&!w)\ eff 

(mnf  (rtf  I2  +  n2w2  )%6  (g+ zf 


EQ.  2.12 
EQ.  2.13 


With  the  elements  as  given  in  the  above  equations,  the  lumped  parameter  model  can  be  repre¬ 
sented  as  shown  in  Fig.  2.7.  Each  series  branch  represents  a  pair  of  spring  and  damping  forces  for 
a  given  value  of  m,  n.  The  number  of  such  branches  depends  on  the  required  accuracy.  The  trunca¬ 
tion  error  can  be  restricted  to  under  1%  with  the  use  of  six  branches. 


*4 - 

Fig.  2.7:  The  spring-damper  lumped  parameter  model 


This  model  has  been  implemented  in  Verilog-A  [8]  for  use  with  the  node-based  MEMS  simu¬ 
lation  tool  developed  at  Carnegie  Mellon  University  [9,  10].  It  has  also  been  implemented  in 
MAST  to  be  used  with  a  similar  library  in  SABER  environment  [11].  The  model  accepts  the  dis¬ 
placement  across  the  terminals  of  the  model  as  lumped  parameter  inputs  and  uses  a  voltage  con¬ 
trolled  voltage  source  to  generate  velocity  by  a  differential  operator.  The  sum  of  the  force  through 
all  the  parallel  branches  represents  the  total  damping  force.  The  model  is  parameterized  and  can  be 
used  in  design.  The  input  parameters  are  length  and  width  of  the  plate,  gap  (the  thickness  of  the 
film)  and  the  ambient  pressure.  The  effective  viscosity  is  calculated  using  Knudsen  number  which 
in  turn  is  evaluated  for  a  given  input  ambient  pressure. 
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Some  equivalent  circuit,  lumped  parameter  models  for  squeeze  film  damping  have  been  pro¬ 
posed  in  the  past  [3,  12].  There  are  a  couple  of  advantages  of  representing  the  model  in  terms  of 
spring  and  damper  elements  rather  than  circuit  elements  like  resistors,  capacitors  and  inductors  [3, 
12].  Firstly,  the  mapping  of  the  through  variable  (force)  and  the  across  variable  (position)  from  the 
mechanical  domain  to  the  electrical  domain  is  not  canonical.  So  is  the  case  with  the  elements  used 
in  the  lumped  parameter  model.  For  example,  the  spring  element  can  be  modeled  by  an  inductive 
or  a  capacitive  element  while  a  damping  element  can  be  modeled  by  a  resistive  or  a  conductive 
(both  dissipative)  element.  Secondly,  in  an  integrated  MEMS  simulation  environment,  there  are 
different  domains  (e.g.  electrical,  mechanical  and  magnetic  domains),  each  of  which  has  its  own 
through  and  across  variable.  If  an  electrical  equivalent  circuit  model  is  used  to  model  mechanical 
damping  elements,  there  arises  a  need  to  separate  the  current  which  represents  the  force  of 
mechanical  domain  from  the  current  that  represents  the  electrical  through  variable.  Hence,  ease  of 
use,  the  use  of  spring-damper  model  is  helpful. 

2.5  Non  idealities 

2.5.1  Edge  and  finite-size  effects 

The  force  equations  EQ.  2.4  and  EQ.  2.5  derived  from  the  linearized  Reynold’s  equation  do 
not  consider  the  edge  and  the  finite-size  effects.  The  edge-effects  are  significant  at  tens  of  microns 
scaled  structure  as  shown  in  section  2.3.  To  study  the  accuracy  of  the  spring-damper  lumped 
parameter  model  based  on  EQ.  2.12  and  EQ.  2.13,  we  compare  the  response  of  the  behavioral 
model  with  that  of  finite-element  simulations,  for  varying  squeeze  film  size.  The  error  in  the  peak 
damping  force  predicted  by  the  behavioral  model  as  compared  to  FEM  simulations  for  a  plate 
oscillating  with  an  amplitude  of  0.1  (am  and  frequency  of  1  kHz  is  shown  in  the  plot  in  Fig.  2.8. 
The  degradation  in  the  accuracy  of  the  model  at  plate  sizes  in  the  range  of  few  tens  of  microns  is 
noticeable.  This  is  because  the  edge  and  the  finite-size  effects  become  significant  at  these  MEMS- 
scaled  plate  sizes. 
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g:  2  nm 

frequency:  1  kHz 
amplitude:  0. 1  pm 


0.  0  -I - , - 1 - 1 - 1 - 1 

0  20  40  60  80  100 

Rate  length  (pm) 

Fig  2.8:  Accuracy  degradation  of  the  model  for  small  plate  sizes 

2.5.2  Edge-effect  model 

There  can  be  several  ways  in  which  the  edge  and  the  finite-size  effects  can  be  modeled.  Here 
we  describe  a  couple  of  approaches  and  outline  our  implementation. 

One  technique  to  model  the  edge  effects  is  the  effective  viscosity  of  the  model  to  compensate 
for  the  underestimation  of  the  damping  force.  This  approach  quantitatively  increases  the  estimated 
damping  force,  but  does  not  provide  a  direct  qualitative  intuition  to  model  the  edge-effects.  As 
effective  viscosity  is  already  used  to  model  the  slip  conditions  encountered  in  narrow  gaps,  a 
model  of  i "\efj= f( slip,  edge )  will  make  it  difficult  for  designers  to  understand  which  effect  is  dom¬ 
inant. 

Another  approach  can  be  to  use  an  effective  pressure  above  the  ambient  pressure  to  take  care 
of  the  non-zero  pressure  at  the  plate  edges.  Pressure  is  one  of  the  input  parameters  of  the  model 
that  is  used  in  the  evaluation  of  the  Knudsen  number  and  the  effective  viscosity.  Hence  this  effect 
is  not  distinct  from  the  effective  viscosity  approach. 

There  is  a  need  to  model  the  edge  effects  using  parameters  that  are  decoupled  from  the  other 
parameters.  The  plate  width  and  length  are  such  independent  parameters.  Hence,  increasing  the 
plate  size  to  account  for  the  non-zero  gauge  pressure  at  the  plate  edges  is  one  straightforward 
method  for  modeling  the  edge  effects.  This  solution  is  also  physically  intuitive  because  the  differ¬ 
ential  pressure  becomes  zero  only  at  a  finite  non-zero  distance  from  the  plate  edge.  So  extending 
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gap:  2  jum 

Oscillation  amplitude:  0.1  (am 
Oscillation  frequency:  1  kHz 


the  plate  edge  and  using  trivial  boundary  conditions  at  the  extended  edge  provide  a  way  of  obtain¬ 
ing  more  realistic  pressure  distribution  across  the  plate.  Fig.  2.9  shows  the  pressure  distribution 
and  decay  beyond  the  plate  edges  when  the  top  plate  is  moving  sinusoidally  with  an  amplitude  of 
0.1  (am  and  frequency  of  1kHz. 

The  non-zero  pressure  distribution  along  the  original  plate  edge  is  accounted  for  by  extending 
the  plate  edge  by  an  amount  5 L  and  using  TBC  at  the  extended  edge  [7]. 

The  value  of  5 L  shown  in  the  Fig.  2.9  illustrates  the  physical  intuition  behind  the  solution 
adopted.  The  quantitative  equation  for  the  plate  extension  is  obtained  by  comparing  the  peak 
damping  force  of  the  original  behavioral  model  with  that  from  FEA.  The  parameters  that  can  affect 
8L  are  identified  as  the  gap,  plate  dimensions,  frequency  of  oscillation  and  the  amplitude  of  oscil- 
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Fig.  2.10:  Linear  fit  with  least  sum  of  squared  errors  for  square  plates  oscillating  with 
frequency  of  1kHz  and  an  amplitude  of  0.1  pm 

lation.  Screening  experiments  were  performed  to  find  the  dependency  of  8L  on  these  factors  by 
varying  these  parameters  one  at  a  time.  All  experiments  were  performed  for  square  plates.  The 
plate  size  was  varied  from  5  pm  to  100  pm,  gap  size  from  1  pm  to  5  pm  and  the  frequency  up  to  1 
MHz.  The  amplitude  of  oscillation  was  kept  under  10%  of  the  gap.  For  above  mentioned  plate 
sizes  and  gaps,  the  value  of  5 L  that  gives  the  closest  match  between  the  behavioral  and  finite  ele¬ 
ment  results  was  determined.  8L  has  a  first-order  dependency  on  the  gap  and  a  second-order 
dependency  on  the  plate  size.  Fig.  2.10  shows  the  plot  of  dUg  as  a  function  of  the  plate  size. 
Larger  plates  cause  more  pressure  perturbation  at  the  plate  edge  and  the  pressure  settles  to  ambient 
pressure  further  away  from  the  edge.  Hence,  larger  plates  require  a  higher  value  of  8L. 

The  best  fit  model  for  8L  is  given  by  EQ.  2.14. 

5L(g,  L)  =  g (0.8792  +  0. 1L)  EQ.  2. 14 

where  L  is  the  plate  dimension  in  microns,  g  is  the  gap  and  §L  is  the  extension. 

When  the  value  of  gap  is  decreased,  we  would  expect  the  squeeze  force  to  increase  as  the  con¬ 
tinuum  squeeze  action  is  proportional  to  1/g3'  This  can  lead  to  the  conclusion  that  8L  should 
increase  when  the  gap  is  decreased.  On  the  contrary,  the  equation  shows  that  8L  is  directly  propor¬ 
tional  to  the  gap.  To  understand  this,  we  need  to  consider  the  regime  of  operation.  We  can  consider 
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the  air  to  be  within  the  continuum  limits  when  the  Knudsen  number  is  «  1  (<0.01 ).  The  molecular 
mean  free  path  of  air  at  standard  temperature  and  pressure  is  0.07  pm.  So  for  gaps  <  7  pm,  the  con¬ 
tinuum  limits  are  no  longer  valid.  Slip  conditions  begin  to  occur.  In  such  conditions,  the  edge 
effects  increase  with  increasing  gap  and  hence  8L  is  directly  proportional  to  gap. 

2.5.3  Extended  model 

The  accuracy  of  the  model  can  be  extended  by  modeling  the  edge  and  the  finite  size  effects. 
The  width  w,  and  the  length  l  of  the  squeeze  film  are  increased  by  5L  to  give  the  effective  width 
w^and  effective  length  leff&s  in  EQ.2.15  and  EQ.  2.16. 

weff  =  w  +  8L(g,  w)  EQ-  215 

le/f=l  +  U(g,l)  EQ.  2.16 

The  values  of  the  spring  and  the  damper  elements  in  the  extended  model  is  given  by  substitut¬ 
ing  weff  and  lejj for  w  and  l  respectively  in  EQ.  2. 12  and  EQ.  2.13.  Upon  extending  the  plate  dimen¬ 
sions,  the  improved  accuracy  of  the  model  is  shown  in  the  plot  in  Fig.  2.11.  The  accuracy  is  within 
5%  for  plate  sizes  in  the  range  of  5  pm  to  100  pm. 


Rate  length  (pm) 


Fig.  2.11:  Improved  accuracy  of  the  model  with  extended  edges  of  square  plates 
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2.6  Squeeze  film  simulation  results 

In  this  section,  squeeze  film  simulation  results  using  the  extended  lumped  parameter  model 
described  in  section  2.5  are  discussed.  Simulation  results  for  the  input  of  three  dominant  harmon¬ 
ics  of  a  square  wave  are  presented.  Behavioral  and  finite  element  simulations  are  compared  for 
varying  aspect  ratio  of  the  plate,  amplitude  of  motion  and  frequency  of  oscillation.  A  system  level 
simulation  of  a  resonator  bandpass  filter  having  squeeze  film  damping  as  the  dominant  source  of 
energy  dissipation  is  presented.  Simulation  results  are  validated  experimentally. 

2.6.1  Step  function  response 

The  squeeze  film  model  has  been  derived  under  the  assumption  of  sinusoidal  oscillations  of 
the  top  plate.  To  study  the  accuracy  of  the  model  for  non-sinusoidal  input  displacement,  behavioral 
simulations  with  the  first  three  harmonics  of  a  square  wave  are  performed  and  the  results  com¬ 
pared  to  the  finite  element  simulations.  Fig  2.12  shows  the  displacement  of  the  top  plate  which  is 
the  sum  of  the  three  dominant  harmonics  of  a  square  wave  of  amplitude  0. 1  jam  and  a  frequency  of 
10  kHz.  The  force  on  an  100  |um  square  plate,  with  a  gap  of  2  pm  and  4  pm,  as  a  function  of  time 
is  plotted  in  Fig.  2.13  (a)  and  Fig.  2.13  (b)  respectively.  The  difference  in  peak  force  is  under  1%  in 
the  case  of  the  2  pm  gap  and  15%  in  the  case  of  the  4  pm  gap. 

2.6.2  Aspect  ratio  dependency 

While  deriving  the  finite  size  and  the  edge-effect  model  and  determining  the  accuracy  of  the 
model  in  the  section  2.5,  square  plates  were  considered.  The  effect  of  aspect  ratio  on  squeeze  film 
damping  is  studied  by  keeping  the  area  constant  and  changing  the  length  to  width  ratio  of  the 

squeeze  film.  Damping  force  on  a  plate  of  area  10000  pm2  is  determined  for  aspect  ratios  in  the 
range  of  1-10.  As  the  aspect  ratio  increases,  as  expected,  the  squeeze  action  decreases  because  air 
underneath  the  plate  with  higher  length  to  width  ratio  can  flow  out  without  much  squeeze  action. 
Fig.  2.14  shows  that  the  behavioral  simulations  match  the  finite  element  simulations  to  within  7%. 
This  validates  the  edge  effect  model  for  higher  aspect  ratios,  even  though  the  model  was  derived 
with  square  plates. 
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Fig  2.12:  Displacement  of  the  top  plate  with  three  dominant  harmonics  of  a  square  wave 


Time  (s) 

Fig.  2.13  (a):  Force  on  a  100  pm  square  plate  with  gap  of  2  pm 


2.6.3  Amplitude  of  motion 

The  damping  model  consists  of  non-linear  spring  and  damper  elements.  Thus,  the  damping 
force  is  not  linearly  proportional  to  the  displacement.  However,  for  displacements  that  are  small 
compared  to  the  gap  (<10%),  the  model  behaves  close  to  a  linear  system. 
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length/width 

Fig.  2.14  (a):  Ratio  of  force  on  rectangular  plate  to  that  on  square  plate  for  a  10"8  m2  area 


length/w  idth 

Fig.  2.14  (b):  Error  in  peak  force  when  compared  to  FEM  simulations 

The  peak  damping  force  for  a  10  kHz  sinusoidal  input  of  amplitude  varying  from  40  nm  to 
200  nm  is  shown  in  Fig.  2.15  (a)  -  Fig.  2.15  (d)  for  square  plates  of  size  20  jam  and  100  pm,  and 
gaps  of  2pm,  4  pm.  The  error,  as  compared  to  the  numerical  simulations  is  within  10%  for  all  the 
four  cases.  The  plot  shows  the  linearity  of  the  system  decreases  with  increasing  gap. 

2.6.4  Frequency  dependency 

Frequency  of  the  plate  oscillations  is  an  important  parameter  in  determining  the  squeeze 
force.  When  the  plate  oscillates  with  low  frequency,  the  air  flow  is  nearly  incompressible.  In  such 
cases,  the  damping  force  that  is  proportional  to  the  velocity  dominates  the  spring  force.  The  total 
squeeze  force  is  proportional  to  the  frequency  for  a  given  amplitude  of  motion  and  for  fixed  phys- 
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0.25  — ♦ —  Behavioral  simulation 


plate  and  2  pm  gap 


Fig.  2.15  (b):  Force  and  error  for  a  20  pm 
square  plate  and  4  pm  gap 


Amplitude  (pm) 


Amplitude  (pm) 


Fig.  2.15  (c):  Force  and  error  for  100  pm  square  Fig  2.15  (d):  Force  and  error  for  100  pm  square 
plate  and  2  pm  gap  plate  and  4  pm  gap 
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0.5  -  Damping  force 


Squeeze  number 


Fig  2.16:  Spring  and  damping  forces  per  unit  displacement  as  a  function  of  squeeze  number 


Fig.  2.17:  Peak  force  for  different  frequencies  for  a  100  pm  square  plate 


Fig.  2.18:  Peak  force  for  a  20  pm  square  plate 


ical  dimensions  of  the  squeeze  film.  At  high  frequencies,  the  flow  invoked  becomes  compressible. 
As  a  result,  the  inertial  effects  of  the  fluid  contribute  to  the  damping.  The  relative  contribution  of 
the  spring  and  the  damping  terms  is  plotted  as  a  function  of  squeeze  number  (EQ.  2.7)  in  Fig.  2.16. 
Cut-off  frequency  is  defined  as  the  frequency  at  which  the  spring  and  the  damping  terms  contrib¬ 
ute  equally.  At  frequencies  comparable  to  and  higher  than  the  cut-off  frequency,  the  inertial  effects 
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are  significant.  As  a  result,  the  gauge  pressure  is  not  negligible  compared  to  the  ambient  pressure 
contradicting  one  of  the  assumptions  in  the  squeeze  film  model  [1].  Hence  the  error  at  the  high  fre¬ 
quencies  can  be  significant.  Fig.  2.17  plots  the  peak  force  as  a  function  of  frequency  for  the  case  of 
a  100  pm  square  plate  with  gaps  of  2  pm  and  4  pm.  For  the  gap  of  2  pm,  cut-off  frequency  for  an 
100  pm  square  plate  is  670  kHz.  From  FEM  simulations,  the  peak  gauge  pressure  at  1  MHz  is  as 
high  as  12  kPa.  This  is  about  12%  of  the  ambient  pressure  (0.1  MPa).  The  assumption  of  small 
gauge  pressure  fails  resulting  in  the  high  error  noticeable  at  frequencies  close  to  1  MHz.  The  error 
in  peak  damping  force  for  a  gap  of  2  pm  is  about  50%  while  it  is  29%  for  a  gap  of  4pm.  This  can 
be  attributed  to  decrease  in  the  squeeze  number  with  the  increase  in  gap.  Fig.  2.18  plots  the  peak 
force  for  a  20  pm  square  plate.  Accuracy  is  within  10%  in  this  case  as  cut-off  frequency  is  16  MHz 
for  g  =2  pm. 
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3.  Lateral  viscous  damping 


3.1  Introduction 

Viscous  lateral  damping,  also  known  as  slide  film  damping,  occurs  when  two  parallel  plates 
are  in  relative  tangential  motion.  The  relative  motion  of  the  parallel  structures  invokes  a  fluid  flow 
that  exerts  a  force  opposing  the  relative  motion.  The  force  exerted  by  this  mechanism  depends  on 
the  velocity  gradient  of  the  fluid  near  the  plane  of  the  oscillating  structure,  the  viscosity  of  the 
medium  and  the  area  of  the  oscillating  structure.  The  velocity  gradient  of  the  fluid  in  turn  depends 
on  the  velocity  of  the  oscillating  plate,  and  its  distance  from  other  parallel  structures  that  are  in  rel¬ 
ative  motion.  Hence  a  slide  film  between  two  parallel  plates  is  characterized  by  the  length  and  the 
width  of  the  plates,  the  distance  between  them,  their  relative  velocity  and  the  viscosity  of  the  film. 


Viscous  damping  is  a  dominant  source  of  energy  dissipation  in  laterally-driven  microstruc¬ 
tures.  A  lateral  resonator  is  a  common  example.  The  proofmass  oscillating  parallel  to  a  substrate 
experiences  viscous  drag  by  fluid  both  on  the  top  and  the  bottom  surfaces.  In  lateral  comb  drives, 
lateral  damping  can  be  a  dominant  source  of  damping.  There  are  four  different  lateral  damping 
sources  in  a  lateral  comb  drive.  Between  the  stator  and  the  rotor  fingers,  there  are  two  slide  films 
associated  with  each  rotor  finger  (one  on  each  side  of  the  rotor  finger)  as  in  Fig.  3.1  (a).  These 
films  are  characterized  by  the  overlap  length  of  the  stator  and  the  rotor  fingers,  the  thickness  of  the 
fingers  and  the  equilibrium  distance  between  the  fingers.  There  is  additional  damping  from  the 
fluid  above  and  below  the  rotor  finger.  This  damping  source  acts  on  the  entire  rotor  area  as  shown 
in  Fig.  3.1  (b).  The  bottom  film  is  characterized  by  the  distance  from  the  substrate  and  the  top  film 
can  be  infinitely  thick  if  exposed  to  atmosphere. 


H  Slide  film  if  fj-jff  J 

8  Rotor  finger 

■  Rotor  truss  fijiiiM  □ 

□  Stator 

Substrate 


(a)  Top  view  showing  the  slide  films  (b)  Cross  section  showing  the  slide  films 

between  the  stator  and  rotor  fingers  on  the  top  and  bottom  of  the  rotor 

Fig.  3.1:  Slide  films  in  a  lateral  comb  drive 


23 


3.2  Slide  film  model 


3.2.1  Model  for  continuum  conditions 


The  oscillations  of  a  planar  structure  in  a  gaseous  medium  invoke  flow  in  the  surrounding 
medium.  This  section  describes  the  theoretical  aspects  and  the  way  the  invoked  flow  is  modeled. 
Consider  a  plate  in  the  X-Y  plane  oscillating  in  the  X-direction.  The  dynamics  of  the  gas  flow 


invoked  can  be  modeled  by  the  one-dimensional  diffusion  equation  (EQ.  3.1)  [13,  14]. 


d  v(z)  dzv(z) 


dt 


E- 


dz‘ 


EQ.  3.1 


where  p  is  the  kinematic  viscosity  of  the  gas  and  v(z)  is  the  gas  velocity  component  in  the  x  direc¬ 
tion.  The  kinematic  viscosity  is  the  ratio  r(/p  where  r|  is  the  viscosity  coefficient  or  the  dynamic 
viscosity  of  the  gas  and  p  is  the  density  of  the  gas.. 


Fig.  3.2:  Slide  film  between  a  planar  oscillating  structure  and  the  substrate 


For  a  fully  established  flow  without  any  transient  consideration,  the  flow  can  be  considered  to 
Couette  flow  with  a  linear  velocity  profile  as  in  Fig.  3.2  (a).  There  are  no  inertial  effects  in  this 
steady  state  (established)  flow.  If  the  continuum  limits  are  considered,  the  velocity  of  the  gas  at  the 
plate  surface  is  equal  to  the  velocity  of  the  plate.  Under  these  conditions,  the  damping  coefficient 
is  given  by  [13,  14].  : 

EQ.  3.2 

9 

where  A  is  the  area  of  the  moving  plate  and  g  is  the  thickness  of  the  slide  film. 

Couette  damping  model  does  not  take  inertial  effects  into  consideration  and  hence  can  be 
inaccurate  at  high  frequencies.  The  frequency  dependent  Stokes  flow  model  takes  into  account  the 
inertial  effects.  For  a  plate  velocity  of  v0cos(cot),  the  velocity  profile  in  Stokes  flow  decays  expo- 
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nentially  as  shown  in  Fig.  3.2  (b)  and  the  phase  shift  varies  linearly  as  we  move  away  from  the 
plate.  Under  continuum  assumptions,  the  gas  velocity  profile  and  the  damping  are  given  in  EQ.  3.3 
and  EQ.  3.4. 

vfz)=  v  sinh(qg-qz)  EQ.  3.3 

0  sinh(qg) 

E,  =  r\A qcoth(qg)  EQ- 

where  q  is  the  complex  frequency  variable  given  by  (/rn/p)  ,  ©  is  the  oscillation  frequency  and  p 
is  the  kinetic  viscosity  and  £,  is  the  complex  damping  coefficient.  A  complex  damping  coefficient 
means  that  there  is  a  damping  force  (the  real  part)  and  there  is  a  spring  force  in  phase  with  the  dis¬ 
placement  (the  imaginary  part).  As  the  frequency  tends  to  zero,  the  Stokes  damping  model 
approaches  the  Couette  damping  model.  The  amplitude  of  the  velocity  profile  decays  by  a  factor  of 
e  over  a  distance  equal  to  the  penetration  depth  5  given  by  (oo/(2p))1/2. 

3.2.2  Gas  rarefaction  effects 

Since  the  critical  film  dimension,  the  gap  between  the  two  plates,  is  in  the  range  of  microns, 
the  mean  free  path  is  comparable  to  the  gap.  This  causes  the  gas  rarefaction  effect,  which  can  be 
modeled  with  an  effective  viscosity.  EQ.  3.5  gives  the  value  of  effective  kinematic  viscosity,  peff 
and  effective  dynamic  viscosity,  r|eff  [14]. 


o  2  4  6  8  10 

Gap  (pm) 


Fig.  3.3:  Effective  kinematic  viscosity  as  a  function  of  gap 
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where  Kn  is  the  Knudsen  number  which  is  the  ratio  of  the  mean  free  path  of  the  gas  molecules  to 
the  distance  between  the  plates.  Fig.  3.3  plots  the  effective  kinematic  viscosity  of  air  at  standard 
pressure  and  temperature,  as  a  function  of  the  plate  distance. 


The  gas  velocity  at  the  plate  surface  is  not  the  same  as  the  plate  velocity  due  to  gas  rarefac¬ 
tion.  Slip  conditions  occur  and  the  velocities  at  the  plate  and  substrate  surfaces,  which  are  the 
boundary  conditions  for  EQ.  3.1,  are  given  in  EQ.  3.6  and  EQ.  3.7  [14]. 

*V<2>-  EQ.  3.6 

z=0 

EQ.  3.7 


v(0)  =  v0+'k- 


dz 


v(g)  =  -X 


dv(z) 


dz 


z=g 


where  vq  is  the  plate  velocity,  X  is  the  molecular  mean  free  path  of  air  and  g  is  the  gap  between  the 
plate  and  the  substrate. 


3.2.3  Solution  for  the  velocity  profile  with  gas  rarefaction 


For  sinusoidal  oscillations  of  the  plate,  the  solution  to  EQ.  3.1  with  boundary  conditions  as  in 


EQ.  3.6  and  EQ.  3.7  is  given  as  [14]: 

VM  =  _V  sinh(c!9  -  qz) + qXcosh(qg  -  qz) 

1  0  (1  +  q2X2 )sinh(qg)  +  2qXcosh(qg) 


EQ.  3.8 


From  the  velocity  profile,  the  damping  coefficient  £,  can  be  extracted  as  in  EQ.  3.8. 


|=Tleff>l  dv(z) 


dz 


.  cosh(qg)  +  qXsinh(qg) 

z=0  eff  (1  +  q2^2  )sinh(qg)  +  2qXcosh(qg) 


EQ.  3.9 


where  r|  is  the  viscosity  coefficient  (also  called  dynamic  viscosity),  A  is  the  area  of  the  plate  sur¬ 
face.  Gas  rarefaction  effects  are  taken  into  account  by  using  peff  and  r)eff  for  ju,  and  r)  respectively 
as  in  equation  3.5. 


As  the  frequency  of  oscillation  tends  to  zero,  the  damping  force  tends  to  that  of  Couette 
damping  force.  For  very  high  frequencies,  this  model  tends  to  Couette  damping  model  with  a  gap 
equal  to  the  molecular  mean  free  path  of  air  (A.). 
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Fig.  3.4:  Lateral  damping  force  on  a  sinusoidally  oscillating  plate 

3.3  Finite  element  simulation  of  slide  films 


FEM  simulations  to  evaluate  the  lateral  damping  force  on  a  plate  oscillating  laterally  above 
an  infinite  substrate  were  performed.  To  verify  the  finite  element  simulations  being  performed,  the 
force  on  a  square  plate  of  size  1mm  on  a  side,  oscillating  with  an  amplitude  of  0.1  pm  and  a  fre¬ 
quency  of  10  kHz  was  determined.  Fig  3.4  shows  the  damping  force  with  time  for  a  gap  of  2  pm. 
The  damping  coefficient  can  be  extracted  from  the  peak  force.  Since  the  gap  is  small  compared  to 
the  penetration  depth  at  10  kHz,  Couette  model  with  linear  air  velocity  profile  between  the  plate 
and  the  substrate  can  be  used  to  evaluate  the  damping  coefficient  as  in  EQ.  3.10  [14]. 

B  =  X^~  EQ.  3.10 


9 


where  g  is  the  gap  between  the  plates,  A  is  the  area  of  the  oscillating  plate  and  T|ej^is  the  effective 
dynamic  viscosity  of  air. 


Fig.  3.5:  Trivial  boundary  conditions  in  finite  element  simulation  of  slide  film  damping 
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The  values  of  damping  coefficient  obtained  from  the  FEA  method  and  the  analytical  method 
are  9.5  pNs/m  and  9.84  pNs/m  respectively.  They  match  to  within  4%  which  verifies  and  quanti¬ 
fies  the  accuracy  of  the  finite  simulations. 

As  in  the  case  of  squeeze  film,  the  type  of  boundary  conditions  used  affect  the  accuracy  of  the 
finite  simulations.  The  use  of  trivial  and  non-trivial  boundary  conditions  in  the  simulations  of  lat¬ 
eral  damping  and  the  advantages  of  each  approach  are  discussed  next. 


3.3.1  Trivial  boundary  conditions 

In  trivial  boundary  conditions  (TBC),  the  slide  film  is  represented  by  a  3-D  structured  block 
enclosed  between  two  wall  surfaces  and  the  edges  of  the  wall  surfaces  are  fixed  at  ambient  pres¬ 
sure.  The  representation  does  not  consider  the  finite  non-zero  thickness  of  the  plate.  To  exploit  the 
symmetry  conditions  for  efficient  simulation,  one  half  of  the  plate  is  used  as  shown  in  Fig.  3.5. 


3.3.2  Non-trivial  boundary  conditions 

The  use  of  non-trivial  boundary  conditions  (NTBC)  provide  the  solver  with  the  freedom  of 
evaluating  the  pressure  and  hence  the  air  flow  around  the  oscillating  structure.  The  oscillating 
plate  is  surrounded  by  blocks  of  air  on  all  directions  and  the  ambient  pressure  boundary  conditions 


i  Substrate 
1  Slide  Film 

I  Oscillating  Plate 


sJU 

Fig.  3.6:  Use  of  non-trivial  boundary  conditions  in  slide  film  damping 
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are  set  at  the  extended  boundaries  to  model  more  realistic  conditions.  Fig.  3.6  shows  the  system 
used  for  the  finite  element  simulations  with  NTBC 

3.3.3  Trivial  boundary  conditions  vs.  non-trivial  boundary  conditions 

The  following  table  provides  the  lateral  damping  force  predicted  by  the  use  of  NTBC  and 
TBC  when  a  plate  of  varying  size  oscillates  with  an  amplitude  of  0.1  pm  and  frequency  of  10  kHz, 
2  pm  above  the  substrate.  It  can  be  observed  from  the  table  that  the  difference  in  force  predicted  by 
TBC,  NTBC  and  the  Couette  damping  model  decreases  as  the  plate  size  increases.  The  dominance 

of  the  edge  and  the  finite-size  effects  for  sizes  less  than  100  pm  on  a  side  has  been  experimentally 
observed  by  Zhang  et  al.  [15]. 


Plate  Size 

Theoretical 

TBC 

NTBC 

(Couette  damping) 

20um 

0.0597  pN/m2 

0.065  pN/m2 

0.0727  pN/m' 

lOOum 

0.0597  pN/m2 

0.0597  pN/m2 

0.064  pN/m2 

1mm 

0.0597  pN/m2 

0.0597  pN/m2 

0.060  pN/m2 

As  in  the  case  of  squeeze  film  damping,  the  force  is  underestimated  when  TBC  are  used.  This 
is  mitigated  by  the  use  of  NTBC,  but  at  the  expense  of  additional  computer  resources.  The  use  of 
NTBC  results  in  the  increase  in  the  number  of  grid  cells  from  100  to  4000  as  a  result  of  which,  a 
typical  transient  simulation  time  increases  from  4  minutes  in  TBC  to  45  minutes  in  NTBC  on  a 
360  MHz  Sun  Ultra-Sparc  machine.  The  memory  requirements  of  the  simulation  is  5  MB  in  TBC 
and  45  MB  in  NTBC.  Non-trivial  boundary  conditions  are  used  for  all  the  simulations  used  to 
characterize  behavioral  model. 

3.4  Lumped  parameter  model 

The  numerator  and  the  denominator  of  the  complex  damping  coefficient  (EQ.  3.9)  can  be  rep¬ 
resented  by  a  Taylor  series  in  variable  q  about  the  dc  value  of  q  (q= 0)  as  in  EQ.  3.11.  The  ratio  has 
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only  even  powers  of  q  in  both  the  numerator  and  the  denominator  which  makes  it  a  ratio  of  com¬ 


plex  power  series  in  angular  frequency. 


5>y 

;~n 


/= 0 


even 


EQ.  3.11 


where  are  the  even  coefficients  in  the  numerator  and  c/(-  are  those  in  the  dominator  and  are  func¬ 
tions  of  the  gap  (g)  and  the  mean  free  path  (X)  of  air  molecules.  The  equations  for  the  coefficients 
are  given  in  the  appendix. 


The  condition  for  the  convergence  of  each  series  can  be  set  by  selecting  the  ratio  of  the  two 
successive  terms  of  the  real  and  imaginary  parts  in  each  series  to  be  less  than  unity.  Truncation 
error  can  be  controlled  by  changing  the  number  of  terms  in  the  series  and  the  ratio  of  successive 
terms.  The  minimum  ratio  of  the  absolute  values  of  successive  terms  in  the  real  and  the  imaginary 

parts  of  the  numerator  and  the  denominator  of  EQ.  3.11  is  g4m2/24[i2.  The  ratio  g4(S)2/24\i2  is  kept 
under  0. 1  to  truncate  the  series  after  first  three  terms.  As  an  upper  bound  on  the  operation  fre¬ 
quency  for  a  given  g  and  p.  For  example  for  g  =  2  pm  and  p  =  17.9  pNms/kg,  co  <  6.9  x  106  s’1. 


The  transfer  function  described  is  implemented  in  Verilog-A  as  a  Laplace  transform.  The  one 
dimensional  slide  film  model  has  two  pins  indicating  the  displacement  of  the  oscillating  plate  and 
the  substrate.  This  is  converted  into  relative  velocity  using  a  time  derivative  operator.  This  relative 
velocity  is  used  to  compute  the  damping  in  the  frequency  domain  and  inverse  Laplace  transform 
back  to  damping  force  between  the  plate  and  the  substrate.  The  input  parameters  of  the  model 
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Displacement 


Length=  2  pm  to  100pm 


Fig.  3.8:  Ratio  of  peak  damping  force  for  varying  length 

include  the  length  and  the  width  of  the  plate,  the  gap  and  the  ambient  pressure.  Fig.  3.7  shows  the 
circuit  implementation  of  the  model. 


3.4.1  Edge  effects  in  the  lumped  parameter  model 

The  damping  force  given  by  EQ.  3.9  does  not  take  into  account  the  edge  and  the  finite  size 
effects  which  are  dominant  at  micron-scale  devices.  This  observation  was  made  experimentally  by 
Zhang,  et  al.  and  it  was  suggested  that  the  damping  force  be  multiplied  by  an  empirical  correction 
factor  of  about  2  to  3  [15].  Since  the  motion  of  the  plate  is  in  the  plane  of  the  structure,  the  edge 
effects  can  show  different  dependencies  on  the  width  and  the  length.  To  study  this,  we  keep  one 
dimension  constant  while  varying  the  other.  Let  the  dimension  of  the  edge  along  the  direction  of 
motion  be  l  (length)  and  that  of  the  perpendicular  edge  be  w  (width).  The  length  and  the  width  are 
varied  one  at  a  time  while  keeping  the  other  fixed  at  100  pm  to  quantify  the  difference  in  damping 
force  obtained  from  behavioral  and  finite  elements  simulations. 

Fig.3.8  plots  the  ratio  of  peak  forces  from  lumped  parameter  and  FEM  simulations  for  lengths 
varying  from  2  pm  to  100  pm.  For  1=2  pm,  the  force  estimated  by  behavioral  simulation  using  the 
model  represented  by  EQ.  3.9  is  less  by  a  factor  of  2.3  than  the  force  from  FEM  simulations,  given 
the  edge  effects  contribution  due  to  finite  width  is  negligible.  This  is  in  agreement  with  the  empir¬ 
ical  factor  used  by  Zhang  et  al.  Peak  force  is  a  sufficient  metric  at  the  frequencies  in  the  range  of 
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Displacement 


Lengths  100pm 


Fig.  3.9:  Ratio  of  peak  damping  force  for  varying  width 

10’s  of  kHz  as  the  damping  fore  dominates  over  the  spring  force  and  the  associated  phase  differ¬ 
ence  is  less  than  1°. 


Fig.  3.9  shows  the  ratio  of  peak  forces  for  varying  width.  The  ratio  of  the  forces  from  behav¬ 
ioral  and  FEM  simulations  is  about  1 .33  for  a  2  pm  width.  This  shows  that  the  edge  effects  depend 
on  the  aspect  ratio  as  well  as  the  orientation  of  the  plate. 

The  contribution  of  the  edge  effects  is  established  qualitatively.  For  quantitative  analysis  and 
derivation  of  the  edge  effect  model,  further  investigation  of  the  dependencies  on  gap  and  the  cou¬ 
pled  effect  of  the  length  and  width  is  required  and  is  left  as  future  work.  Further  simulations  of  the 
lateral  damping  model  will  be  presented  with  a  fixed  plate  extension  of  4  pm  on  each  side.  This 
translates  into  a  multiplication  factor  of  about  1.08  for  a  100  pm  square  plate  which  is  close  to  the 
ratio  shown  in  Fig.  3.8  and  3.9. 


3.4.2  Step  function  response 

The  damping  model  used  in  the  derivation  of  the  lumped  parameter  model  was  derived  for 
sinusoidal  input.  To  study  the  response  to  non-sinusoidal  input,  a  sum  of  the  first  three  dominant 
harmonics  of  a  square  wave  of  amplitude  0.1  pm  was  used  as  input.  Fig.  3.10  are  plots  of  the  lat¬ 
eral  displacement  and  force  on  a  100  pm  square,  oscillating  plate.  The  plate  dimensions  were 
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extended  by  4  |im  for  which  the  damping  force  gets  multiplied  by  a  factor  of  1.08,  which  is  close 
to  the  edge  effect  ratio  in  Fig.  3.8  and  Fig.  3.9. 

3.4.3  Frequency  dependency 

The  lateral  damping  force  is  highly  frequency  dependent.  At  low  frequencies,  the  fluid  near 
the  plate  surface  can  respond  to  the  plate  oscillations.  In  this  case,  the  phase  difference  between  the 
force  and  the  velocity  is  zero.  Inertial  effects  are  negligible  and  hence  there  is  no  spring  compo¬ 
nent.  As  the  frequency  increases,  inertial  effects  come  into  play  and  the  force  on  the  oscillating 
plate  has  both  a  damping  and  a  spring  component.  The  total  force  can  be  analyzed  by  comparing 
the  amplitude  and  phase  of  the  force  from  the  behavioral  and  the  FEM  simulations.  Fig.  3.11  (a) 
shows  the  plot  of  magnitude  as  a  function  of  frequency  on  a  100  pm  plate.  It  can  be  seen  that  the 
peak  force  is  almost  linear  with  frequency,  at  least  until  1MHz,  after  which  the  amplitude  increases 
more  rapidly  with  increasing  frequency.  The  phase  difference  between  the  force  and  the  velocity  is 
plotted  in  Fig.  3.11  (b).  The  phase  is  close  to  zero  for  frequencies  up  to  10  kHz.  Phase  increases 
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linearly  with  frequency.  The  phase  response  of  the  behavioral  model  matches  to  within  2°  to  that 
of  finite  element  model  for  frequencies  up  to  1  MHz. 


1000  10000  100000  1000000 

Frequency  (Hz) 

Fig.  3.11  (a):  Force  amplitude  on  a  100  pm  square  plate 
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4.  System  level  simulation  and  experimental  verification 

This  chapter  provides  the  system-level  simulation  of  a  bandpass  filter  and  a  folded-flexure 
microresonator.  The  simulation  is  done  using  the  NODAS  cell  library  developed  at  Carnegie  Mel¬ 
lon  University  [9,  10].  The  bandpass  filter  has  squeeze  film  as  dominant  damping  source  and  the 
lateral  resonator  has  slide  film  damping  as  the  dominant  source. 

4.1  CMOS  bandpass  resonator  filter 

To  emphasize  the  use  of  the  behavioral  squeeze  film  model  in  capturing  the  dynamics  of 
damping  accurately,  results  of  a  system  level  simulation  of  a  bandpass  resonator  filter  are  pre¬ 
sented  and  experimentally  validated. 

Fig.  4.1  shows  the  schematic  diagram  of  a  550  kHz  CMOS  micromechanical  bandpass  filter 
[16].  The  filter  is  designed  using  the  NODAS  parameterized  cell  library.  The  filter  system  is  com¬ 
posed  of  three  resonators  coupled  by  “O”  springs  and  is  integrated  with  the  associated  electrical 
interface  circuitry.  Each  resonator  is  designed  for  a  550  kHz  resonance  frequency.  The  system  has 
three  resonant  peaks  scattered  around  the  natural  frequency  of  550  kHz.  The  filter  consists  of  three 
pairs  of  differential,  electrostatic  comb  drives  used  for  actuation  and  sensing.  The  differential 
comb  drives  have  associated  squeeze  film  damping  which  is  the  dominant  source  of  energy  dissi¬ 
pation  in  the  comb  drives  as  well  as  in  the  resonator  filter.  The  squeeze  film  between  the  stator  and 
the  rotor  fingers  has  a  length  of  37  pm,  a  width  of  4.8  pm  and  a  thickness  of  1.8  pm.  The  number 
of  such  squeeze  films  in  each  comb  drive  is  twice  the  number  of  rotor  fingers  (one  on  each  side  of 
the  rotor  finger).  The  lateral  damping  model  used  in  the  simulation  is  a  Couette  flow  with  the  gap 
assumed  equal  to  the  penetration  depth  at  550  kHz. 

The  output  voltage  obtained  at  the  sensing  comb  drive  is  buffered  by  the  electronics.  Fig  4.2 
shows  the  output  voltage  measured  at  the  output  of  the  interface  circuitry.  With  the  use  of  the 
damping  model  with  edge  effect  enhancements,  the  output  voltage  at  the  resonant  frequencies 
matches  the  experimental  measurements  to  within  0.6  dB,  which  is  equivalent  to  7%  error  in  the 
gain.  The  results  obtained  with  the  use  of  the  behavioral  model  without  any  edge  extensions  differ 
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Driving  Coupling  Coupling  Sensing 

resonator  resonator  spring  resonator 

(b):  Bandpass  resonator  filter  schematic 
Fig.  4. 1 :  NOD  AS  schematic  of  the  bandpass  resonator  filter 
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Fig  4.2:  AC  transfer  characteristics  of  the  bandpass  filter 
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Folded  flexure 


Anchor  points 


Fig.  4.3.  Layout  view  of  the  lateral  folded-flexure  comb-drive  microresonator 

by  8  dB,  a  gain  error  of  over  150%.The  lateral  damping  associated  with  the  filter  is  smaller  than 
the  squeeze  damping  by  over  an  order  of  magnitude. 

4.2  Lateral  microresonator 

In  a  lateral  microresonator,  the  lateral  slide  film  damping  can  be  a  dominant  energy  dissipa¬ 
tion  mechanism.  Fig.  4.2  shows  the  layout  view  of  a  folded-flexure  lateral  resonator  [17].  The  sys¬ 
tem  can  be  modeled  by  a  mass-damper-spring  system  characterized  by  a  resonant  frequency  and 
quality  factor. 

The  slide  films  associated  with  the  all  the  moving  components  like  the  comb  drive,  the  beams 
and  the  shuttle  mass  have  associated  slide  film  damping  contributing  to  energy  dissipation  in  the 
system.  Different  slide  films  in  the  comb  drive  include  the  films  between  the  stator  and  the  rotor 
fingers  and  between  the  rotor  and  the  substrate.  A  detailed  description  of  these  is  provided  in  sec¬ 
tion  3.1.  At  frequencies  less  than  100  kHz,  the  damping  force  on  the  bottom  surface  of  the  micro¬ 
structure  dominates  as  the  distance  between  the  substrate  and  the  structure  (2  pm)  is  much  less 
than  the  penetration  depth.  At  frequencies  close  to  or  greater  than  100  kHz,  the  penetration  depth 
of  air  becomes  less  than  10  pm.  Under  such  conditions,  the  damping  force  on  the  top  surface  con¬ 
tributes  significantly  to  the  overall  damping. 
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Fig.  4.4:  Frequency  response  of  a  7.4  kHz  microresonator 


As  mentioned  in  chapter  3,  a  constant  edge  effect  model  has  been  used  in  which  all  the 
dimensions  are  extended  by  4  (4m.  Fig.  4.4  shows  the  plot  of  resonator  displacement  as  a  function 
of  frequency.  The  damping  associated  with  the  top  surface  is  not  included  as  the  penetration  depth 
at  10  kHz  frequency  is  about  26  pm.  This  makes  the  damping  associated  with  the  top  surface  over 
an  order  of  magnitude  smaller  than  the  damping  on  the  bottom  surface.  This  can  contribute  signif¬ 
icantly  for  frequencies  close  to  1  MHz.  Fig.  4.4  shows  the  frequency  response  of  a  7.4  kHz 
microresonator  obtained  from  NOD  AS  simulation.  Table  4.1  shows  simulated  and  experimental 
values  of  resonant  frequency  and  quality  factor  for  this  resonator.  Simulation  results  are  in  good 
agreement  with  the  experimental  results.  The  simulated  resonant  frequency  matches  to  within 
1.3%  and  the  simulated  Q  factor  matches  to  within  8.5%  of  the  experimental  values.  Table  4.1  also 
shows  the  results  of  microresonators  designed  for  30  kHz  and  100  kHz  resonant  frequency.  It  can 
be  seen  that  the  error  in  the  Q  value  increases  with  increasing  frequency.  This  can  be  attributed 
mainly  to  two  reasons.  The  constant  edge  extension  model  is  an  approximate  model  for  the  edge 
effects.  It’s  accuracy  needs  to  be  verified  for  frequencies  close  to  100  kHz.  The  damping  on  the  top 
surface  (exposed  to  atmosphere)  that  is  not  taken  care  of  in  the  schematic  contributes  to  the  overall 
damping. 

Table  4.1:  Resonant  frequency  (fr)  and  quality  factor  (Q)  for  different  microresonators 


design  frequency 

fr  (experiment 

Q  (experiment) 

fr  (NOD  AS) 

Q  (NODAS) 

10  kHz 

7.4  kHz 

8.3 

7.31  kHz 

9 

30  kHz 

22.7  kHz 

25.6 

21.9  khz 

31.6 

100  kHz 

74.1  kHz 

80 

72.6 

130 
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5.  Conclusions  and  future  work 

This  thesis  describes  behavioral  models  for  squeeze  film  and  slide  film  damping  in  MEMS 
devices.  The  models  are  parameterized  and  can  be  used  in  mixed  domain  system  level  simulation 
of  MEMS.  It  is  shown  that  the  edge-effects  contribute  significantly  to  the  overall  damping  for 
plate  sizes  <50  pm  on  a  side.  Extension  of  the  plate  edges,  as  one  of  the  possible  ways  to  model  the 
edge-effects,  is  discussed.  The  edge  effect  model  is  derived  for  the  case  of  squeeze  film  for  gaps  in 
the  range  of  1  pm  to  5  pm,  plate  sizes  in  the  range  of  5  pm  to  100  pm.  For  slide  films,  it  is 
shown  that  the  edge  effects  are  dependent  on  both  the  aspect  ratio  and  the  orientation  of  the  slide 
film.  Gas  rarefaction  effects  have  been  modeled  using  an  effective  viscosity  model. 

To  show  the  utility  of  the  models  in  modeling  the  damping  dynamics,  results  of  system  level 
simulations  have  been  presented  and  experimentally  validated.  With  these  models  the  designer 
will  be  able  to  model  the  squeeze  film  and  the  lateral  (slide  film)  damping  with  adequate  accuracy. 

There  are  several  directions  for  future  work.  The  presence  of  holes  in  the  microstructure  sig¬ 
nificantly  alters  the  damping  force  on  the  oscillating  structure.  The  effect  of  the  hole  size  and  the 
hole  density  on  the  damping  force  needs  to  be  studied.  A  constant  edge  extension  model  has  been 
used  to  model  the  edge-effects  in  lateral  damping.  The  gap  and  plate-size  dependent  edge  effect 
model  derivation  can  extend  the  accuracy  of  the  edge-effect  model  to  wider  range  of  input  param¬ 
eters.  The  effect  of  the  thickness  of  the  oscillating  plate  on  the  squeeze  film  and  lateral  damping 
needs  to  be  assessed.  The  study  of  the  effect  of  the  presence  of  other  static  or  moving  structures 
near  a  squeeze  film  or  a  slide  film  can  be  useful  in  understanding  the  damping  in  case  of  complex 
geometries. 


39 


References 


[1]  J.  J.  Blech,  “On  Isothermal  Squeeze  Films”,  Journal  of  Lubrication  Technology,  V.  105,  1983, 
pp. 615-620. 

[2]  J.  B.  Starr,  “Squeeze  Film  Damping  in  solid  State  Accelerometer”,  Tech  Digest,  IEEE  Solid 
State  Sensor  and  Actuator  Workshop,  Hilton  Head  Island,  SC,  USA,  June  1990,  pp.44-47. 

[3]  T.  Veijola,  H.  Kuisma,  J.  Lahdenpera,  and  T.  Ryhenen,  “Equivalent-Circuit  Model  of  the 
Squeezed  Gas  Film  in  a  Silicon  Accelerator”,  Sensors  and  Actuators  A,  vol.  48,  1995,  pp.239- 
248. 

[4]  CFD-VIEW,  CFD-GEOM,  CFD-GUI  User  Manuals,  CFDRC,  Huntsville,  AL. 

[5]  M.  Torowski,  Z.  Chen,  A.  prezekwas,  “High  Fidelity  and  Behavioral  Simulation  of  Air  Damp¬ 
ing  in  MEMS”,  Proc.  MSM  ‘99,  pp.  241-244. 

[6]  E.  S.  Kim,  Y.  H.  Cho,  and  M.  U.  Kim,  “Effect  of  holes  and  edges  on  the  squeeze  film  damping 
of  perforated  mico-mechancial  structures”,  IEEE  Int.  Conf.  MEMS,  1999,  pp.  206-301. 

[7]  S.  Vemuri,  G.  K.  Fedder,  T.  Mukherjee,  “Low-Order  Squeeze  Film  Model  of  Simulation  of 
MEMS  Devices”,  Proc.  MSM,  2000. 

[8]  Verilog-A  reference  Manual. 

[9]  J.  E.  Vandemeer,  “Nodal  Design  of  Actuators  and  Sensors  (NODAS)”,  M.S.  Thesis, 
Carnegie  Mellon  University,  1998. 

[10]  G  K.  Fedder,  Q.  Jing,  “A  Hierarchical  Circuit-level  Design  Methodology  for  Micormechani- 
cal  Systems”,  Trans.  Circuits  and  Systems  II,  TCAS  1999. 

[11]  MAST  Reference  Manual,  Release  4.2,  Analogy  Inc.,  Beaverton,  OR,  1997. 

[12]  T.  Veijola,  T.  Ryhanen,  H.  Kuisma,  and  J.  Lahdenpera,  "Circuit  Simulation  Model  of  Gas 
Damping  in  Microstructures  With  Non-Trivial  Geometries,"  Transducers ’95  -  Eurosensors  IX, 
Stockholm,  June  1995. 

[13]  Y.  H.  Cho,  B.  M.  Kwak,  A.  P.  Pisano,  and  R.  T.  Howe,  “Slide  film  damping  in  laterally  driven 
microstructures”,  Sensors  and  Actuators  A,  vol.  40,  1994,  pp.  31-39. 

[14]  T.  Veijola,  “Compact  Damping  Models  for  Lateral  Strucutres  Including  Gas  Rarefaction 
Effects”,  Proc.  MSM ,  2000. 

[15]  X.  Zhang,  and  W.  C.  Tang,  “Viscous  Air  Damping  in  laterally  Driven  Microstructures”,  Sen¬ 
sors  and  Materials,  vol.  7,  No.  6,  1995,  pp.  415-430. 

[16]  Q.  Jing,  H.  Huo,  T.  Mukherjee,  L.  R.  Carley,  and  G  K.  Fedder,  “CMOS  Micromechanical 
Bandpass  Filter  Design  using  a  hierarchical  MEMS  Circuit  Library”,  IEEE  lnt.  Conf.  MEMS, 
2000. 

[17]  T.  Mukherjee,  S.  Iyer,  G.  K.  Fedder,  “Optimization-based  synthesis  of  microresonators”,  Sen¬ 
sors  and  Actuators  A  70,  1998,  pp.  118-127. 


40 


Appendix 


The  equation  for  slide  film  damping  coefficient  including  the  gas  rarefaction  effects  is: 

e  _  Tleffd  dv(z)  =  cosh( gg)  +  gXsinhf qg) 

v  dz  z=Q  eff  (1  +  q  X  )sinh(qg)  +  ZqKcosh(qg) 

The  numerator  and  the  denominator  are  odd  functions  in  the  complex  frequency  variable  q 
and  hence  have  only  odd  powers  in  their  taylor  series  expansion  in  q  about  q=0. 

Let  the  numerator  be  N(q)  and  the  dominator  be  D(q).  The  taylor  series  expansion  of  N(q)/ 
(rleffA) is: 

q+q3(g2/2+gX)+q5(g4/24+g3A/6)+q7(g6/720+g5A/120)+q9(g8/40320+g7A/5040)+q11(g1°/ 

3628800+g9A/362880)+0[q13]. 

Similarly  the  denominator  D(q)  is  given  by: 

q(g+2A)+q3(g3/6+g2A.+gA,2)+q5(g5/120+g4y\./12+g3X2/6)+q7(g7/5040+g6A/360+g5A2/ 

1 20)+q9(g9/3 628 80+gsA/20 1 60+g7A2/5040)+q  J1(gn/3991 6800+g  joa/18  14400+g9A2/ 
362880)+O[q13]. 

The  ratio  of  the  two  series  has  only  even  powers  of  q  in  the  numerator  and  denominator  as  in: 

oo 

2"*' 

:=0 

f-  — 

S  oo 

i=0 

even 

The  Coefficients  n;  and  dj  are  as  in  the  above  mentioned  series. 
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